%% author:sousanone
%% email:nj2000515@foxmail.com
%% date:2023/4/30
%% 文献参考：Extended coherent modulation imaging for single-shot object retrieval free from illumination artifacts



clc,clear all;

addpath 'D:\matlab2022\bin\cmi'

% addpath('D:\mat\bin\ptychography-master\utils')
%% FZA generation
lamda = 632.8e-6;                                        % 入射波长，单位mm米、, 632.8nm nm-um-mm-m
k0=2*pi/lamda;                                        
Nx =600; 
range1=2000;
focallength=30;                                         %焦距 
step1=range1/(Nx-1);                                     % 步距
yp1=-range1/2:step1:range1/2;
xp1=-range1/2:step1:range1/2;
[x1,y1]=meshgrid(xp1,yp1);                               % 划分网格
position1  = x1 + 1i * y1;                               % 某点的位置
lensletr = abs(position1);                               % 某点的半径
lensang  = angle(position1);
Nc=floor(Nx/2)+1;
TFZP= exp(-1i*pi*lensletr.^2/(lamda*focallength));
a=sign(real(TFZP));
lens1 =sign(1+sign(real(TFZP)));   



%% CMI simulation
lamda = 632.8e-6;
pix = 5.5e-3;
z1 = 20;
z2 = 10;
z3 = 15;
% % mask1=zeros(Nx,Nx);
% % mask1(Nx/2-200:Nx/2+200,Nx/2-200:Nx/2+200)=1;

% % mask1=circle_mask(Nx,400,Nx/2,Nx/2);


lens2=lens1;
% % % 这啥也不是
% cmi_prior = lens2;%.*exp(1i*lens2/max(max(lens2)));

% % 相位型
cmi_prior = exp(1i*lens2/max(max(lens2)));

% % 振幅形
% cmi_prior = lens2/max(max(lens2)).*exp(1i*mask1);


obj=imread('cameraman.tif');
obj=imresize(obj,[Nx,Nx]);
obj1=double(obj);
phase=imread('mla.jpg');
phase=double(rgb2gray(phase));
phase=imresize(phase,[Nx,Nx]);
phase=phase;

% % obj2=double(obj1).*exp(1i*phase/max(max(phase)));


%%
sup=circle_mask(Nx,400,300,300);
probe= f_Gaussian_beams(1e-3,0,lamda,2e-3,2e-3,600,600);
pphase=imread('Lena.tif');

pphase=double(imresize(pphase,[Nx,Nx]));
probe=double(probe).*exp(1i*pphase/max(max(pphase)));
%% probe乘以支撑区域约束，前向传播到物体处
ew1=double(sup).*probe;
iw1=Propagate(ew1,'fresnel',pix,lamda,z1);
%% 物体与照明光结合形成出射波ew2
obj2=phase.*exp(1i*obj1/max(max(obj1)));
ew2=iw1.*obj2;
%% 出射波ew2前向传播到mask部分，然后乘以mask形成出射波ew3
iw2=Propagate(ew2,'fresnel',pix,lamda,z2);
ew3=iw2.*double(cmi_prior);
%% 出射波前向传播到CCD处，得到入射波iw3
iw3=Propagate(ew3,'fresnel',pix,lamda,z3);
raw=abs(iw3);







iter = 15;
obj_g = double(ones(Nx,Nx));
guessprobe=double(im2bw(abs(probe)));


for i = 1:iter

    %% 猜测的probe乘以支撑区域约束，前向传播到物体处
    new1=guessprobe.*sup;
    niw1=Propagate(ew1,'fresnel',pix,lamda,z1);

    %% 猜测的物体与照明光结合形成出射波new2
    new2=niw1.*obj_g;
    niw2=Propagate(new2,'fresnel',pix,lamda,z2);

    %% 出射波new2前向传播得到iw2到达mask部分，然后乘以mask形成出射波new3
    new3=iw2.*double(cmi_prior);
    if i>=2
        new3=2*new3-temp;%%论文中提到的更新方式。
    end

    %% 出射波前向传播到CCD处，得到入射波niw3
    niw3=Propagate(new3,'fresnel',pix,lamda,z3);
    
    %% 进行幅值替换
    dew1= raw.*exp(1i.*angle(niw3));

    %% 逆传播到MASK处
    diw2 = Propagate(dew1,'fresnel',pix,lamda,-z3);
    if i>=2
        diw2=diw2-new3+temp;%%论文中提到的更新方式。
    end
    %% 更新一下经过mask的出射波dew2
%     dew2 = diw2 + conj(double(cmi_prior))./(max(max(abs(double(cmi_prior)).^2))).*(diw2-new3);
    %动量项更新方式
    b1=0.5;
    dew2 = diw2 + conj(double(cmi_prior))./((1-b1).*abs(double(cmi_prior)).^2+b1*(max(max(abs(double(cmi_prior)).^2)))).*(diw2-new3);
% %    


    %% 反向传播到物平面
    diw3 = Propagate(dew2,'fresnel',pix,lamda,-z2);

    %% 更新一下物光波和照明波前
%     obj_g= obj_g+conj(niw1)./max(max(abs(niw1).^2)).*(diw3-new2);
%     guessprobe1=niw1+conj(obj_g)./max(max(abs(obj_g).^2)).*(diw3-new2);

    


    b2=01;
    b3=0.04;
    obj_g= obj_g+conj(niw1)./((1-b2).*abs(niw1).^2+b2*max(max(abs(niw1).^2))).*(diw3-new2);
    guessprobe1=niw1+conj(obj_g)./((1-b3).*abs(obj_g).^2+b3*max(max(abs(obj_g).^2))).*(diw3-new2);


    guessprobe=Propagate(guessprobe1,'fresnel',pix,lamda,-z1);
    guessprobe=guessprobe.*sup;
    temp=new3;




    if mod(i,10)==0
        %figure,imshow(abs(ew1),[])
    end
end
% subplot(221)
% imshow(abs(obj_g),[])
% title('物体振幅')
% subplot(222)
% imshow(angle(obj_g),[])
% title('物体相位')
% subplot(223)
% imshow(abs(guessprobe),[])
% title('probe振幅')
% subplot(224)
% imshow(angle(guessprobe),[])
% title('probe相位')


subplot(221)
imshow(abs(obj_g),[])
title('物体振幅')
subplot(222)
imshow(angle(obj_g),[])
title('物体相位')
subplot(223)
imshow(abs(guessprobe),[])
title('probe振幅')
subplot(224)
imshow(angle(guessprobe),[])
title('probe相位')


